function [gk] = grobalStiffnessMatrix(coords,ne,nodes,nelx,nely,e,nu,phi)
%Construct global stiffness matrix

lambda = e*nu/((1+nu)*(1-2*nu));
mu = 0.5*e/(1+nu);

%ke = reshape(elementStiffnessMatrix(e,nu),64,1);
phib = max(0.001,min(1.0,phi.^3));

idx = reshape([1 2 3 4 5 6 7 8]'*ones(1,8),64,1);
jdx = reshape(ones(8,1)*[1 2 3 4 5 6 7 8] ,64,1);

is = zeros(ne*64,1);
js = is;
ss = is;
l = 1;
for i = 1:ne
    nd = nodes(i,:);
    vtx = coords(nd,:);
    edofs = [nd(1)*2-1;nd(1)*2;nd(2)*2-1;nd(2)*2;nd(3)*2-1;nd(3)*2;nd(4)*2-1;nd(4)*2];
    is(l:l+63) = edofs(idx);
    js(l:l+63) = edofs(jdx);
    
    ke = reshape(stima4(vtx,lambda,mu),64,1);
    
 %   area = 0.5*(det([1 1 1; vtx(1:3,:)']) + det([1 1 1; vtx(2:4,:)']))

    ss(l:l+63) = phib(i)*ke;
    
    l = l + 64;
end
gk = sparse(is,js,ss,2*(nelx+1)*(nely+1),2*(nelx+1)*(nely+1));

end

